home *** CD-ROM | disk | FTP | other *** search
/ WINMX Assorted Textfiles / Ebooks.tar / Text - Mathematics - Numerical Mathematics and Computing (F).zip / bvp2.f < prev    next >
Text File  |  2002-06-11  |  2KB  |  68 lines
  1. C
  2. C PAGE 421-420: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
  3. C
  4. C FILE: BVP2.FOR
  5. C
  6. C BOUNDARY VALUE PROBLEM SOLVED BY SHOOTING METHOD (RK4SYS,XPSYS)
  7. C
  8.       DIMENSION  X(5),X2(99),X4(99)   
  9.       DATA  NP1/99/, TA/1.0/, TB/2.0/ 
  10.       DATA  ALPHA/1.09737491/, BETA/8.63749661/ 
  11.       DATA  X/1.0,1.09737491,0.0,1.09737491,1.0/
  12.       H = (TB - TA)/REAL(NP1) 
  13.       T = TA      
  14.       DO 2 I = 1,NP1
  15.         CALL RK4SYS(5,T,X,H,1,0)      
  16.         X2(I) = X(2)
  17.         X4(I) = X(4)
  18.         T = TA + REAL(I)*H  
  19.    2  CONTINUE
  20.       P = (BETA - X4(NP1))/(X2(NP1) - X4(NP1))  
  21.       Q = 1.0 - P 
  22.       DO 3 I = 1,NP1
  23.         X2(I) = P*X2(I) + Q*X4(I)     
  24.    3  CONTINUE
  25.       ERROR = EXP(TA) - 3.0*COS(TA) - ALPHA     
  26.       PRINT 5,TA,ALPHA,ERROR
  27.       DO 4 I = 9,NP1,9      
  28.         T = TA + REAL(I)*H  
  29.         ERROR = EXP(T) - 3.0*COS(T) - X2(I)     
  30.         PRINT 5,T,X2(I),ERROR 
  31.    4  CONTINUE
  32.    5  FORMAT(5X,F10.5,2(5X,E22.14))   
  33.       STOP
  34.       END 
  35.         
  36.       SUBROUTINE XPSYS(X,F) 
  37.       DIMENSION  X(5),F(5)  
  38.       F(1) = 1.0  
  39.       F(2) = X(3) 
  40.       F(3) = EXP(X(1)) - 3.0*SIN(X(1)) + X(3) - X(2)      
  41.       F(4) = X(5) 
  42.       F(5) = EXP(X(1)) - 3.0*SIN(X(1)) + X(5) - X(4)      
  43.       RETURN      
  44.       END 
  45.         
  46.       SUBROUTINE RK4SYS(N,T,X,H,NSTEP,J)
  47.       DIMENSION  X(10),Y(10),F1(10),F2(10),F3(10),F4(10)  
  48.       H2 = 0.5*H  
  49.       START = T   
  50.       DO 6 K = 1,NSTEP      
  51.       CALL XPSYS(X,F1)      
  52.       DO 2 I = 1,N
  53.    2  Y(I) = X(I) + H2*F1(I)
  54.       CALL XPSYS(Y,F2)      
  55.       DO 3 I = 1,N
  56.    3  Y(I) = X(I) + H2*F2(I)
  57.       CALL XPSYS(Y,F3)      
  58.       DO 4 I = 1,N
  59.    4  Y(I) = X(I) + H*F3(I) 
  60.       CALL XPSYS(Y,F4)      
  61.       DO 5 I = 1,N
  62.    5  X(I) = X(I) + H*(F1(I) + 2.0*(F2(I) + F3(I)) + F4(I))/6.0     
  63.       T = START + REAL(K)*H 
  64.    6  IF(K .EQ. (K/J)*J)  PRINT 7,T,(X(I),I = 1,N)
  65.    7  FORMAT(5X,F10.5,4(5X,E22.14))   
  66.       RETURN      
  67.       END 
  68.